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I. INTRODUCTION 


A successful engineering design is usually a compromise between the capacity of 
the elements in the design and the actual operating conditions. These elements as 
well as the operating conditions are usually subject to such random effects that they 
become random variables. Hence, the properties of the elements involved and the 
external conditions are not known with complete certainty. Two extreme results of a 
classical engineering design are those related to under- or over-design which led to 
failure of the design or to a large waste of resources. It is essential to guarantee that 
no catastrophic limiting states occur during the operation or lifetime of a design. 
Thus, the extreme values are the ones affecting the failure of any given system. 
Consequently, the problem of the distribution of extreme value is of great interest. 

The pioneering work of E. J. Gumbel which was published in his book “Statistics 
of Extremes” in 1958 had a large impact on the engineering community. Since that 
time many new advances and developments have taken place in the area of extreme 
value theory. These new developments are in engineering applications as well as in 
areas of probability, reliability, stochastic processes, multivariate structures, and 
statistical decision theory. The engineering applications involve a broad range of 
topics such as breaking strength of material [3], extreme river stages [3,4], wind shear 
[1], wind energy [8], aircraft gust loads [12], etc. However, with few exceptions, most 
of these applications have been restricted to the univariate cases. 

An engineering application, which is of particular interest to us, is the “ascent 
structural wind loads” for the National Space Transportation System (NSTS). The 
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response of the space vehicle to loads induced by the winds, in particular, the 
fluctuation in the structural loads due to wind profile variability, has been and 
continues to be the subject for intensive investigations. An aerospace vehicle must 
withstand the design ascent structural load or be within any flight constraints to 
safely reach orbit. For the NSTS an assessment for flight readiness is made for each 
mission. This assessment for ascent loads and performance parameters is stated in 
terms of launch probability, which is derived from analyses of certain probability 
models. 

The univariate extreme forms are the Gumbel type, which is unbounded, the 
Weibull (of largest extremes) and the Frechet, which have lower bounds. Smith et al. 
[13,14] used the bivariate Gumbel “m~case” distribution in the analyses of extreme 
loads for the Space Transportation System (STS). However, some of the load 
variables had a physical lower bound of zero, and hence application of the bivariate 
Gumbel was not strictly appropriate. Therefore, the need for practical models for the 
bivariate Weibull and Frechet was clearly identified for this specific application. This 
is the motivation for this report. It is envisioned that in many other engineering 
fields, in which the univariate models of extremes are used, applications of the 
bivariate models will evolve after dissemination of the tools provided in this report. 

Extreme value distributions were obtained as limiting distributions of greatest 
(or smallest) values in random samples of increasing size. Since 

min X: = -max(-Xj), minima results can be obtained from the corresponding 

1 < i < n 1 < i < n 

maxima, and vice versa. 

Let (X A , Yj), i = 1, 2, •*, n, be a sequence of i.i.d. random pairs each having the 
same continuous joint cumulative distribution 

F(x, y) = P{X i < x ,Y i <y}. 
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We consider the joint distribution of X max := mao^Xj, X 2 , • • •, X n ), and 
Y max := max(Y 1 , Y 2 , •••, Y n ). 

To obtain a nondegenerate limiting distribution, it is necessary to reduce the 
actual greatest value by applying a linear transformation whose coefficients depend on 
the sample size. That is, we find linear transformations 

X( n ) &n X max + b n , ^>0 , and Y( n ) C n Y max “b * C n >0. 
such that each of and Y^ n j has a limiting distribution which is one of three types 
of extreme value distributions (see section II below). The limiting joint distribution 
° f X (nb Y (n) is a bivariate extreme value distribution. The joint cumulative 
distribution function of X max and Y max is [F(x,y)] n . Denoting the bivariate extreme 
value cumulative distribution function by F^^y), we then have 

F (oo)( x .y) = lim FKx + b n, c n y + d n )] n , (1) 

For notational convenience we shall write F(x,y) for the bivariate extreme 
distribution, if no ambiguity is likely to arise. 

Equation (1) is sometimes referred to as the “stability relation” or “stability 
postulate.” It is a natural extension of a similar equation obtained by Frechet, and 
also used by Fisher and Tippett, in the univariate case. 
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II. UNIVARIATE EXTREME VALUE DISTRIBUTIONS 

The univariate extreme value distributions are among the three following 
families; 

(1) Gumbel: P{X < x} = exp [-e“ x ], 

(2) Frechet: P{X < x} = { ^ ^ ^<° 

(3) Weibul!: P(X<x} = { «xp [ - ( - x) ^ x ^ 0 „ 

where X — -- a - is a reduced variate, x (bold x) is the original variate, /i is the 
location parameter, a is the dispersion parameter, and 6 > 0 is a parameter. The 
corresponding distributions of ( -X) are also called extreme value distributions. 

Of these three families of distributions, Gumbel is by far the most commonly 
referred to in discussions of extreme value distributions, see [3]. 

The Weibull distribution above (sometimes called “reversed” Weibull) is an 
asymptotic distribution of smallest extremes. It is used in the modeling of failure (or 
survival), and those applications connected with reliability analysis. 

The mirror image of the Weibull distribution, the distribution of (-X), is called 
Fisher-Tippet type III distribution of largest extremes, 

(3)' P{ - X > - x} = 1 -P{-X<-x} = l - exp ( -x)° , x > 0. 

The distribution in (3)' has been used to model maximum temperatures, maximum 
wind speeds [7], and maximum earthquake magnitudes. 

The Frechet distribution is an asymptotic distribution of largest extremes 
derived by Frechet under the condition that the initial variates be nonnegative. It 
has been widely used as a model for floods and maximum rainfall [2,7]. 
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II.l. REMARK: We note that the Frechet and Weibull distributions, with parameter 
9 , can be transformed into Gumbel distribution (with and ^=0) by the 

transformations Z=lnX, and Z= -ln( -X), respectively. Therefore, we can obtain 
the Frechet or the Weibull distributions, with parameter 9 —\ > 0, from Gumbel type 
with /z=0 and <r=^ , by using the transformations z=e x or z=-e' x , respectively. 
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III. STRUCTURE OF BIVARIATE EXTREME VALUE DISTRIBUTIONS 

Gumbel and Mustafi [5] have described two general forms for the bivariate 
extreme value distributions in terms of the marginal distributions; 

(1.) The a -case: 

F(x,y) = F,(x) F 2 (y) exp[-a(j^ + j^y)' 1 ], (1) 

where 0 <a< 1, X and Y are reduced variates. The parameter a indicates the 
association between the two extremes, and a=0 corresponds to the case where X 
and Y are independent. 

(2.) The m -case: 

F(x,y) = exp[- {In F^x))™ + ( -In F 2 (y)) m }®], (2) 

where 1 < m < oo, X and Y are reduced variates. The parameter m indicates the 
association between the two extremes, and m=T gives the independence case. 

We shall consider bivariate extreme value distributions when the margins F 1 (x) 
and F 2 (y) belong to the same family of univariate extreme value distributions. 

A. Gumbel Margins 

J. Tiago de Oliveira [15,16] proved that a bivariate extreme value distribution 
with Gumbel margins can be defined by a cumulative distribution function of the 
form 

F(x,y) = exp[ - (e‘ x + e' y ) g(y - x)], (3) 

where X and Y are reduced variates. 

In this case a density function exists if the “dependence” function g satisfies the 
following conditions: 

<*) few = fe. « (t) = '• 

(ii) Tt K 1 + O §(01 < 0, 
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(iii) £ [(1 + e ‘) g(t)] > 0, and 

(iv) (1 + e _t ) g ;, (t) + (l-e"*) g*(t) > 0. 


A.l. Remarks: (i) To obtain the a -case equation (1) from equation (3) above, we 

take 

g(t) = gi(t) = 1 - f sech 2 (§), 0 < a < 1. 

(ii) To obtain the m-case equation (2) from equation (3), we take 


(4) 


g(t) = g a (t) = “It 0 - . l<m<oo. 


e +1 


We will now verify these remarks. 

(i) The a - case: 

F(x,y) = exp[ - (e" x + e‘ y ) (1 - - * (y . x) )] 


4cn 


= exp[ - (e" x + e" y ) (1 - f (— 2 —j^) 2 )] 

e 2 -j- e ^ 

= exp[ - (e' x + e " y ) + a + \J x) ~ 1 

(e 2 -f e 2 


(5) 


(y + x) 


= exp[-(e' x + e' y ) + a 


. y-x _ (y~x) J 

2 T e 2 J 


= exp[ - (e’ x + e" y ) + a e y . 


Therefore, 


F(x,y;a) = exp[-( e - x + e‘ y ) + a (e x + e y n, 0 < a < 1. 


( 6 ) 


(ii) The m - case: 

F(x,y) = exp[ - (e" x + e' y ) 


(e m(y-x) + 


.y- x 


+ i 


= exp[ - (e x 


+ e' y ) ( s 


-mx + e ‘ my 
e' my 
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Therefore, 


F(x,y;m) — exp[-(e~ mx + e' my ) m ], 1 < m < oo. 


( 7 ) 


A.2. Rotation of the Axes by 45° 

We may want to take advantage of the symmetry around the diagonal y=x by 
rotating the coordinates 45°. 


x = x , cos 45° - y' sin 45° - ^ (x 1 - y'). 
y = x' sin 45° + y' cos 45° = ^ (x*+ y ; ). 


(i) The a - case: 


F(x',y’) = exp[-(e'¥ (x '- , ' ) + + a 

^2 J , "fl I ^ 

'] 


e f( x '+ y, )Vn 


= exp [ - (e 2 * (e 2 y -f e 2 * j j 1 - ^ 


T2 ^2 / 

( •'*?•’ 7 


V2 J 


= exp[- 2 e 2 ch(^ y') (1 - § sech 2 (4^ y'))], 0 < a < 1. 


(8) 


Thus, 


V2 , 

X 


F(x',y') = exp[- 2e 2 x ch(^y') g^y')], 


where g x (t) = 1 - f sech 2 (M , 0 < a < 1, is the dependence function in A.l. 


( 9 ) 


(ii) Them -case: 

F(x',y') = exp^e^*'^ e^ x ' +y,) f\ 




^2 / nc/2 ; _nW2 ,1 
= exp[-e 2 ^e 2 y + e 2 y ^ n ] 


T2 I 


= exp[ - e 2 


x’ / 


^2y ' 

W~i 


+1 


m-g-y 




1 mfy' ) 


= exp(-e*>^ + l) ] 
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A.3. The Joint Density functions 

(i) The a - case: 

f(x,y;a) = F(x,y;a) e' (x+y) [l-a(e 2x + e 2y )(e x + e y )' 2 + 2 a e 2(x+y) (e x + e y )' 3 

+ a 2 e 2(x+y) (e x + e y )‘ 4 ] (11) 

where F(x,y;a) is the distribution function given by (6), and 0 < a < 1. 

(ii) Them -case: 

f(x,y;m) = F(x,y;m) e ' m(x+y) (e mx + e ‘ my )~ 2+ R f [m-1 + (e mx + e‘ my )% (12) 
where F(x,y;m) is the distribution function given by (7), and 1 < m < oo. 

A.4. The Bivariate Mode 

(i) The a -case: 

The value (x , y ) which maximizes the density f(x,y;a), given by (11), is called 
the modal location. This value lies on the line y = x (or y ; = 0), and hence to obtain 
it we write the density in terms of x' and y', the 45° -rotated axes. From (8), we 
have, 

x + y = >J2 x' , e 2x + e 2y = 2 e^ 2 x ( e — y e ^~ y ) = 2 e^ x ch(>|2 y') , and 
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£x' 

e x -f e y = 2 e 2 ch(^- y'). 


Therefore (11) becomes 

f(x',y';a) = F(x',y';a) - a 2 y '> j. 


+ 


(2 e 2 X ch(f y')) 2 (2 e 2 X ch(^ y')) 3 


a 2 e 2 ^ 2 *' 1 

+ 2 * 2 , r- J ' 

(2 e 2 ch(^y')) 4 


f(x , ,y / ;a) = F(x',y';a) e'^fl -§ ch(^ y') sech 2 (^y') + f e 2 x sech 3 (^ y') 


+ sech 4 (^P y')] - 


Thus, 


V2 , 


f(x',0;a) = F(x',0;a) e'^ x '[l -§ + fe + f e 2 X 


_V2 , 

where F(x',0;a) = exp [-2(1 - f )e ] 
Therefore, 
f(x',0;a) — exp 


-2( 1 - | )e'^ x '] ((l-f) 2 ^' + fe^'), 


(13) 


and 


^2 , 


< 2 , 


^-f(x',0;a) = exp [-2( 1 - | )e^ ] (1 ((1 -f) 2 ^' + f ) 


^2 , 


^ , 


+ exp [ — 2(1 — f)e 2>c ]( — >[ 2 (i — |) 2 e'^ x/ - ^ e"" 2 ”* ) . 
Setting ^-f(x',0;a) = 0 , we obtain 

(1 -|) e-^*' ((1 -f) 2 e^*' + f e'? x ') - (1 -f) 2 e-^ x ' - § e^' = 0 


(1 — f) ((i-|) 2 + 


^2 , 


^ , 


^2 0 -Af2x' x‘e'2 1 ). (! -|) 2 e - 2 * -| = 0 


(l-f)V^ -(l-f)(l-f)e- 2 
Therefore, 


^ , 


= 0 


-fx' (1 - |)(1 - §) ± J(X - f) 2 (l -f) 2 -4(1- 1) 3 ( - |) 


2(1 - 1) 3 
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Therefore, modal location is y'= 0, and 



The corresponding modal value is obtained by substituting (14) into (13): 
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(1^2 i mij2 / 


= e 


( e 2 + e 2 y ) 


m\f2 i 

T~* 


i 

= 2 e 2 x ch(^y'). 


Therefore, (12) becomes 


Af2 / 


f(x',y';m)=F(x',y , ;m)e' m ^ x ^2e 2 x ch(2^ y')) 2+I " [m - 1 


+ 2e? x W^y#]- 


(16) 


Thus, 

f(x',0;m) = F(x',0;m)e" mN ^ x ^2e 


^2 


f(x',0;m) = exp[-2 ,n e 2 x ] 2 2+m ^(m -l)e 2 x + 2 m e x> j 



mT2 / 

(2e 

2 x \ 

^2 

/ 

2 5 

], or 

t 

x ] 


/ 

x ] 

2 - 2 +ri f( 


. . >f2 i 

y + ™ (m-l+ 2 e 2 X ), 


)' 2+m [m - 1+ 2 e 


_£ x / 
2 x 


V2 , j 

~ V 1 


(17) 


Now, differentiating w.r.t. x' and setting the derivative to zero, we obtain 
^-f(x',0;m) = exp[ - 2™ e 2 * ] (-2™ (^)e 2 x ) ((m - 1 ) e 2 x + 2 m e ” ^ 


i ^f 2 J f ^f 2 J i | — / \ 

+ exp[ - 2 m e 2 (m - i) e 2 + 2 rn (-^ T )e “ ^ x ) =0, 


which gives, 


JL -^x'/ 

2 rn e 2 ^(m - i) e 


^2 


? x ' + 2^ e"^ x ') -((m - i) e ^ x + 2^ e -^ x ') 


2^( m -3)e~^ x '+ 2^ e 2 X 


¥ x i 

(m - 1 ) e 2 =0 


2m e ->f2 x ' + 2 m ( m -3) e 2 


^ , 
X 


- (m - 1 ) =0. Thus, 
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^2 f 
2 x 


- 2 m (m - 3) ± \2 m (m - 3) 2 + 4(2) m (m - 1) 

y 

2(2) m 


- (m -3) 


i (m — 3) + 4(m — 1) 


2(2) m 


m > 1. 


It follows that, 


^2 t 
2 x 


(m — 3) -f ^ (m — 3) + 4(m — 1) 

( 2 )m + i 


- (m - 3) +>] (n 

(2) f 


■ l) 2 + 4 


+ 1 


, and 


— ^y X ; = ln[ — (m — 3) +>J (m — l) 2 + 4 ] — +1) ln2. 

Thus the modal location is y' = 0, and 

x'= 42 {(£+») ln2 - ln[ - (m - 3) (m-l) 2 + 4] }. 

Substituting (18) into (17) we obtain the modal value 


(18) 


f(x' ,0;m) = exp[- 


(m — 3) +-^ (m — l) 2 + 4 . / — (m — 3) +, (m — l) 2 + 4 


-!(- 


(( m -!) + (■ 


— (m — 3) +,J (m — l) 2 + 4 


)) 


exp[- 


(m — 3) + ^ (m 


l) 2 + 4 


K- 


(m — 3) +^J (ni — 1 ) 


+ 4 ^(m+1) +>J (m ■ 


l) 2 + 4 


) ( 19 ) 
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B. Frechet Margins 

Using the transformation z = e x or In z = x, (we use Gumbel type with /i=0 
and <r=i see Remark II.l.) in equations (3), (6), (7), (9), (10), (11), and (12), we 
obtain the following (after changing the variables back to x and y): 


B.l. Bivariate Extreme Distributions (with Frechet margins) 

The equation that corresponds to equation (3) is 
F(x,y) = exp[-(x‘ < ’+ y' B ) g((9(ln y-ln x))] 

= exp[-(x' <> + y~ 6 ) g{9 In })], > 0, x > 0, y > 0, (20) 

and F(x,y) = 0, for x < 0, y < 0, 
where X and Y are the original variates (with /z = 0). 

The function g is the dependence function given in Remarks A.l. The choice of 
g, as gj or g 2 , will yield (respectively): 


(i) The a - case: 


F(x,y;a) = |°’ 

|exp[-(x 


if x < 0, y < 0, 
e + y* S ) + a(x fl + y 0 )' 1 ], if x > 0, y > 0, 


where 9 >0 is a parameter and 0 < a < 1. 
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B.2. Rotation of the Axes by 45° 

(i) The a -case: 


F(x',y';a) = j °’ 


if x' < 0, y ; < 0 , 


^2 


(23) 


exp[-2(x')'^ch(f0 In y')(l -f sech 2 (^ In y'))], if x' > 0 , y' > 0 


< 2 , 1 „/\ 


where 9 > 0, and 0 < a < 1. 


(n) The m -case: 
F(x',y';m) = J 0 ’ 


ex 


if x' < 0 , y' < 0 , ^ 


where 9 > 0, and 1 < m < oo. 

Note that the two cases above (equations (23) and (24)) can be written as 

F(x',y')={ 0 ' x' < 0 , y' < 0 , {J5) 

1 exp[ - 2(x') ^ * ch(f 0 In y') g(ln(y')' Bl ') ], if x' > 0 , y' > 0 , 


where 9 > 0, and g is as given in Remarks A.l. 


B.3. The Joint Density Functions 

(i) The a -case: 

f(x,y;a) = F(x,y;a) !>:y) "(1 _ a(x 2 » + y M )(x' + y")' 2 + 2a(xy) M (x ,i +y V 3 

+ a 2 (xy) 2 V + y # )‘*l. ( 26 > 

where F(x,y;a) is the distribution function given by (21), # > 0, 0 < a < 1, x > 0, 
y > 0, and f(x,y,;a) = 0, otherwise. 


(n) The m -case: 

f(x,y;m) = F(x,y;m) (xy ) ~ m6 (x _me + y _n,e ) -2 + nf [m-l + (x _md + y~ m *)™ ] (27) 

where 6 > 0, 1 < m < oo, x > 0, y > 0, and f(x,y;m) = 0, otherwise. 
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B.4. The Bivariate Mode 

(i) The a-case: 

To write the density function (26) in terms of x' and y', the 45°-rotated axes, 
using relations (8), we have xy = i(x' 2 -y' 2 ) , x 20 + y 20 =( ^ ) 20 [(x' -y') 26 

+(x'+y') 2<? ], and x e + y 6 = ( [(x' -y') 0 +(x'+y')® ]. 

Now, the Frechet density function (26) takes the form 

f(x',y';a) = Ffx',y':a)( x * ~ y V * [1 -ail - y') 29 +(x'+y') 29 )((x' -y'f +(x'+y') 9 )' 2 

+ 2a (f)' 3 " «x'-y') 9 +(x'+y') 9 )~ 3 

+ a 2 ((x'-y') 9 + (x' + y')r ], 

= F(x',y';a)2» (x' 2 -,' 2 )- 9 [l-a((x'-y') 29 +(x'+y') 29 )((x'-y')“ +(x'+y')V 

+ 2 ‘-?a(x' 2 -/ 2 ) 29 ((x'-y') 9 +(x'+y') 9 )' 3 

+ a 2 (x' 2 -,' 2 ) 29 ((x'-y') 9 +(x'+y')V], 
where F(x',y';a) is as given (23), 9 > 0, 0 < a < 1 , x' > 0 , y' > 0 , 
and f(x',y';a) = 0, otherwise. 


(28) 


The modal location for the Frechet can be obtained from the corresponding value 
for the Gumbel, equation(14), as follows yp = e y = e° = 1, and 


Xf = >f2 exp{ J [In (2(1 -f) 2 ) - In ((1 -|)+\|l - f + x)]} 


= a|2 ((2(1 - f) 2 )" ((1 - |) + y|l- f + £)'"], 
where 0= i . 


(29) 


rv 

The modal value is, then, obtained by substituting y' F and x F , given above, 
into equation (28). 
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(ii) The m-case: 

As in the a-case, the density (27), can be written as 

f(x',y';m) = F(x',y';m)(x' 2 -y' 2 )- m(? (2)5((x'-yT m<? + (x'+yT 1 "*)' 2 ^ 

x( m-1 + (2)5 ((x'-y')' 1 "* + (x'+y')' 1 "* )® ), (30) 

where F(x / ,y / ;m) is given by (24), 8 > 0, 1 < m < oo, x' > 0 , y' > 0, and 
f(x',y';m) = 0, otherwise. 

The modal location is obtained as in the a-case; yp = 1, and 
Xf = >f2 exp[ \ {(i+1) ln2 - ln(-(m-3) (m-1) 2 + 4 )}] 

= 42 (2)^ +1) (-(m-3) (m-1) 2 + 4 )'*, (31) 

where 6 = ^ . 


C. Wei bull Margins 

For notational convenience we present here bivariate extreme distributions of 
smallest values, the bivariate distributions of largest extremes can be easily obtained 
from it, see equation (3)' of section II and the discussion therein. Using the 
transformation z = -e~ x or ln(-z)=-x, we obtain: 


C.l. Bivariate Extreme Distributions (with Weibull margins) 

F(x,y) = 


exp[ - (( - x) 0 + ( - y) 0 ) g (8 In f)], x < 0, y < 0, 
0, x > 0, y > 0, 


(32) 


where > 0, and X,Y are the original variates (with //= 0). 

The dependence function g is as given in Remarks A.l. 

(i) The a-case: 


F(x,y;a) 

<.-1 


_ j exp[-((-x) e + (-y) e ) + a((-x)’ # + ( -y) )'*], x < 0, y < 0, 


0 , 


x > 0, y > 0, 


( 33 ) 


where 0 =^ >0 is a parameter and 0 < a < 1. 
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(ii) Them -case: 

F , v v ._* = f exp[-((-x) me + (~yr d )h x<0, y<0, 
*Tx,y,mj x > 0, y > 0, 

where 8 =& >0 is a parameter and 1 < m < oo. 


(34) 


C.2. Rotation of the Axes by 45° 

(i) The a - case: 


exp[-2( -x')~ 2 < ch(^0 ln(-y'))(l-f sech 2 (^ ln(-y'))], 




F(x',y';a) = 


0 , 


x' < 0, y' < 0, 
x' > 0, y' > 0, 


(35) 


where 8 > 0, and 0 < a < 1. 


(ii) The m - case: 


F(x',y';m) = Q 


exp[-2(-x')’ 2 S ch(^ ln(-y')) (( * ] ^ + x' < 0, y' < 0 


((-y ) ^ +i) / 36 n 

x' > 0, y' > 0, { > 


where 0 > 0, and 1 < m < oo. 

Again, we notice that both the cases above can be written as 

^ 2 , 


F(x',y') = 


exp[-2(-x')~ e ch(?0 ln(-y')) g(ln( -y')^)], x' < 0, y' < 0, 


0 , 


x' > 0, y' > 0, 


(37) 


where 8 > 0, and g is as given in Remarks A.l. 


C.3. The Joint Density Functions 

(i) The a -case: 

f(x,y;a) = F(x,y;a) (xy)‘ [l-a((-x)-“ + (-y)‘ M ) ((-x)-» + (-y)'V 

+ 2a (xy)' M ((-x)- < ’+ ( -y) ®r 3 + a J (xy)- M (( -x)-»+ (-y)-')' 4 ], (38) 

where F(x,y;a) is the distribution function given by (33), 8 > 0, 0 < a < 1, x < 0, y < 0, 
and f(x,y;a) = 0, otherwise. 
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(ii) Them -case: 

f( x ,y;m) = F(x,y;m) (xy) m0 (( - x) m0 + (-y) m0 )' 2 + i ^ 

x [m-1 + ((-x) m0 + (-y) m0 )^], (39) 

where F(x,y;m) is the distribution function given by (34), 0 > 0, 1 < m < oo, x < 0, 
y < 0, and f(x,y;m) = 0, otherwise. 


C.4. The Bivariate Mode 

The bivariate mode for the Weibull is similarly obtained as that of the Frechet. 
(i) The a-case: 

The Weibull density function (38) takes the form 
f(x',y';a) = F(x',y';a)2 0 (x' 2 -y' 2 )’ 0 [1 -a((x'-y') 20 + (x / +y , ) 2fl )((x'-y')*+ (x'+y') 0 )' 2 
+ 2 1- ^ (- 1) 0 a(x' 2 -y' 2 ) 20 (( x '- y ') 0 + ( x '+ y ') 0 )' 3 
+ a 2 (x' 2 - y /2 ) 20 «x'-y') 0 +(x'+y')V], (40) 

where F(x',y';a) is as given in (35), 0 > 0, 0 < a < 1, x' < 0, y' < 0, and f(x',y';a) = 0, 
otherwise. 


The modal location for the Weibull can be obtained from the corresponding 
value for the Gumbel, equation(14), as follows y(y — -e" y = — e® = — 1, and 

x w = - 'B exp( % [In (2(1 -S) 2 ) - ln((l - }) + \|l - | + t)]1 

= -'12 1(2(1 -f) 2 )'* ((1 -f) + \|l - f + 4 )'], (41) 

where 0 = 

The modal value is, then, obtained by substituting y' w and x(y, given above, 
into equation (40). 
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(ii) The m-case: 

As in the a-case, the density (39), can be written as 
f(x',y';m) = F(x',y';m)(x' 2 - y' 2 )’ 1 "* (2)* -yV 00 + (x'+y')' 1 "*)' 2 ^ 

x (m - 1 + (2)5 (-l) fl ((x / -y , )" mfl + (x'+y'p*)®) , (42) 

where F(x',y';m) is given by (36), 6 > 0, 1 < m < oo, x' < 0, y' < 0, and f(x',y';m) = 0, 
otherwise. 

The modal location is obtained as in the a-case; = — 1, and 
x(v = — 42 exp[ ^ {(m +1) ln2 - ln(- (m-3) + >J(m-l) 2 + 4 )}] 

= - 42 (2)- e ^ +1) ( - (m - 3) + (m - l) 2 + 4 )*, (43) 

where 9 — 

The modal value is, then, obtained by substituting y' w and x' w , given above, into 
equation (42). 
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IV. PROPERTIES OF BIVARIATE EXTREME DISTRIBUTIONS 

Bivariate extreme value distributions have several properties (for example, see 
Gumbel [3], and Tiago de Oliveira [19]); we only discuss a few here. 

(1.) For any bivariate extreme value distribution the following holds 

F {oo) (x,y) > F 1(oo) (x) F 2(oo) (y) , 

where Fj^^x), and F 2 ^ 00 j(y) are the (marginal) univariate extreme distribution 
functions. 

(2-) K F(.)(* ,y) and G^foy) are two bivariate extreme distributions, then so is 
their weighted geometric mean 

( F (0 .,(x, y ))*(G M ( x ,y)) , -\ 0 < A < 1. 

(3.) Frechet bounds. Let F(x,y) be any bivariate distribution function with 
marginals Fj(x) and F 2 (y). Then 

max (0, F 1 (x)+F 2 (y) - 1) < F(x,y) < min (F^x) , F 2 (y)). 

(4.) Each of F(-,«;a) and f(-,-;a), the bivariate extreme distribution and density 
functions, given by equations (6) and (7) respectively, is an increasing function of the 
parameter a , see Gumbel and Mustafi [5]. 

(5.) The Gumbel type distribution, Fj^dx) = exp(-e~ x ), has expected value 

2 

7 (=0.5772- • •), the Euler’s constant, and variance ^L-. 

(6.) Since lim (e' mx + e"^) 1 " = max(e' x , e' y ), it follows that 

m-+oo v 7 \ 7 / ? 

lim F(x,y;m) = min[exp( -e' x ), exp( -e" y )]. 
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(7.) The conditional distributions: We shall use the formula 

F *'Y WY=y) =tey (1 > 

to drive the conditional distributions. 

A. Gumbel margins 
(i) The a -case: We have 

F(x,y;a) = exp[-(e' (ii) * * * * * * * x -f e‘ y ) + a (e x + e y ) _1 ], 0 < a < 1, 

thus 

~ = ex P[-( e " x + e" y ) + a (e x + e y ) _1 ] (e~ y - a e y (e x + e y )‘ 2 ). 

Also, 

37 F y(y) = 37 ( ex P(- e ' y )] = exp(-e' y ) (e' y ). 

Using equation (1) above, we obtain 

F x | Y (x 0 |Y=y) = exp[ -e x ° + a(e x ° + e y ) _1 ] (l-ae 2y (e >0 + e y )' 2 ) 

= ( J - a( TT^) 2 ) exp[ -e _x ° + a(e x ° + e y ) _1 ] (2) 

where 0 < a < 1, X and Y are reduced variates. 


(ii) The m-case: We have 

F(x,y;m) = exp[ - (e" mx + e‘ my )™], 1 < m < oo, 

= exp[-(e‘ mx + e' my )®] (_^(e' mx + e' my )= _1 (-m) e' my ) 
= exp [ - (e‘ mx + e' my )^] (e~ mx + e' 1 ” 5 ')^ _1 e' my . 


Also, 

37 F y (y) = exp[ -e" y ] (e‘ y ). 

Equation (1) yields: 

F x | Y (x 0 |Y=y) = (e mx °+ e' my )™ _1 exp[ - (e mx °+ e' my ) iK - (m - l)y + e' y ], 
where 1 < m < oo and X,Y are reduced variates. 

If we let z = e mx °+ e" my , equation (3) takes the form 

F X I Y ( x ol Y= y) = 1 exp[-x™ -(m-l)y + e' y ]. 


( 3 ) 


( 4 ) 
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Equation (3) can also be written as 


? x i y ( x ol Y =y) = 


i , 

/ " mx 0_L o m y\m ' 1 

| ) exp[ - (e' mXo + e m *)™ + e' y ] 

(m-l)y 


= ( e~ mX °+ e^ )- u - 1 exp[ _ (e' mx o + e -n. y) m + e -y] (5) 

where 1 < m < oo, X and Y are reduced variates. 

Now, using the transformations given in section II, we can obtain the 
corresponding formulas for the Frechet and Weibull margin cases. 


B. Frechet margins 
(i) The a -case: 

F x | Y ( x o|Y=y) = (l-a(^^) 2 )exp[-x 0 <, + a(x* + yV], ( 6 ) 

where O<a<l,0 = ^>O, and x 0 > 0, y > 0 are values of the variates X and Y. 


(ii) The m - case: 

F x , Y (x 0 |Y=y) = ex p[ - ( x o m6 + y' me r + y' 9 ]. ( ? ) 

where 1 < m < oo, 0 = ^ > 0, and x 0 > 0, y > 0 are values of the variates X and Y. 

C. Weibull margins 

(i) The a -case: 

F x | Y ( x ol y =y) = l 1 _ J-e+^y)-* ) 2 ) ex p[ - ( - x °)'+ a (( - x o)'+ (-y)'VJ* (8) 

where O<a<l,0 = ^>O, and x 0 < 0, y < 0 are values of the variates X and Y . 

(ii) Them -case: 

F x|Y (XolY=y)=( ( ~ x ° > ^ ) ( J ) — f" exp[-((-x 0 ) m# + (-yrt*+ (-y)*l, 0) 

where 1 < m < oo, 0 = ^ > 0, and x 0 < 0, y < 0 are values of the variates X and Y. 
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V. PARAMETER ESTIMATION 

There are five parameters for each type (Gumbel, Frechet, and Weibull). These 
are the location and dispersion (shape or scale) parameters fi and a for each of the 
two variates and the association parameter a or m as the case may be. 

A. The location and dispersion parameters 

The marginal reduced variates are 
Zj — (X - /q) and Z 2 = ^ (Y- /x 2 ). 

The estimates i=l,2, are the modal values, and the values of the shape 
parameters <tj, i=l,2, can be calculated from the separate marginal distributions 
using any of the following methods: 

1. The method of maximum likelihood. 

2. The method of moments. 

3. Best linear unbiased estimates (Lieblein’s method). 

4. Best linear invariant estimators. 

The method of maximum likelihood yields a system of equations whose exact 
solutions are difficult to obtain. However, Kimball developed a simpler procedure to 
get approximate solutions. This is described by Kimball in some sections of 
Gumbel’s book, pp. 229-235 of [3]. 

The method of moments, as proposed by Gumbel [3], for the variate X gives 
(and similar corresponding results for Y) 

Zj = X - , where X 1 is a reduced variate. 

Thus, 

*1 7 = X-zq or fi x = X -7 <r u 
where 7 (= 0.5772) is Euler’s constant. 
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Also, 


&ff = var(x) = sj or = $ s x , 


where S* = E . 

i = 1 

If the sample size N is small 7 and are replaced by the corresponding 
population parameters 77 N and cr N , which are tabulated as a function of the sample 
size N. This is referred to as u Gumbel modified moment method.’ 1 In this case, we 
have 


— o 

£1 = X - Sr x £ n and ^ . 

Lieblein proposed to estimate a parameter 9 by 0, a linear combination of order 
statistics in which weights are determined to yield minimum variance unbiased 
estimators. In the 4 th method of best linear invariant estimators, the condition of 
unbiasedness was relaxed from Lieblein’s method. 

Gumbel’s moment method is the most simple of all four methods. 


B. The a and m parameters 

There are three possible methods that can be used to estimate the parameters a 
and m. The first method is based on the contents of the quadrants (the quadrant 
method). The second method is based on relating the parameters to a correlation 
coefficient. The third method (the strip method) is based on the distribution of 

l x -y|- 

We shall discuss these methods briefly, but first we need the following relations. 
F(x,x;a) = exp [-2 e‘ x + § e" x ] 

= exp[-(2-§) e' x ] 

= exp[-e’ (x ‘ ln(2 ^ ) ] 

= F(x-ln( 2 -|)) 

= [F(x)f 1 , 


(1) 
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where F(x) is the first asymptotic (univariate) distribution (F(x)= exp( -e" x )) and 
0<a< 1. 

F(x,x;m) = exp( - 2 m e’ x ) 

= [F(x)]»“ (2) 

where F(x) is the first asymptotic distribution and 1 <m<oo. 

To obtain the median x we set 
F(x) = exp( -e‘ x ) = \ 

=> e" x = -ln(|) or x= -ln(-ln|) 

Thus x = y = 0.3665129. 


B.l. The Quadrant method (see Gumbel and Mustafi [5]) 

(i) The a - parameter 

The population medians of the marginal 
distribution divide the domain of variation of 
x and y into 4 quadrants as shown. The 
number of sample points, Q t and Q 3 , 
contained in quadrants 1 and 3 are obtained. 

The “critical content” is defined by 

Q + Q 

3 = — , where N is the total number 
of sample points. 

Note that with increasing a, i.e. increasing association, the distribution tends to 
concentrate along the diagonal. Therefore, the maximum likelihood estimate of a is 
the value a such that, (from equation (1)). 

F(x,x;a) = (e' e = Q 13 when x=x . 


(reduced vanaifi') 

y! 


a. 




Q. 


<*,a> 


x. 


(reduced vaciftts) 

Q, 
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Since 

F(x,x;a) = (±) 2 ^ 

= 1 ( I)' 1 

= \ (2) f , 

it follows that 
4 Q m = 2* 
or 

In 4 + In Qi ? 3 = | In 2 

2 In 4 + 2 In Q, ~ 

1 >_ ^ — A 

In 2 

A , 2 ln ^1,3 _ . 

4 + TTT"- a - 

Because of the fact that 0 < a < 1, we see that 
F(x,x;0) = (i) 2 2 = \ , and 

F (x,x;l) = (i)5 = 0.35355. 

Therefore, if Qj 3 lies outside the interval [0.25, 0.35355], the quadrant method cannot 
be used to estimate the parameter a. 


( 3 ) 


(ii) The m - parameter 

From equation (2), we have (a+x=x) 

F(x,x;m) = (i) 2 "*. 

Since 1 < m < oo, we see that 
F(x,x;l) = = i and 

F(x,x;oo) = (I) 2 ° = \ . 

Therefore, if Q 13 is outside the interval [^,|] the quadrant method cannot be used 


to estimate m. Otherwise, as in the a -case, 
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Qi , 3 = ar 

_ -L 

ln Q 13 = -2™ In 2 

ln( - In Q 13 ) = -| In 2 + In In 2 


ln(-ln Qj 3 ) - in(ln 2) 
_ 


(4) 


Note that if Q 1 3 is such that 0.25 < Q 13 < 0.35355, both a and m can be estimated 
by the quadrant method. 


B.2. fi) Estimating the parameter a using the grade correlation coefficient 

The grade correlation coefficient was introduced by Hoeffding [6]. Let t be the 
difference sign correlation coefficient, k' be the (Spearman) rank correlation coefficient, 
and k be the grade correlation coefficient. If Z 1? Z 2 , Z n , n > 2, is a random sample 
from a bivariate population F (see I\onijn[9]), where X i =(X i ,Y i ), then 

(n-l)(n + l) _ J_ 


3n* 


-L Esgn(X i -X j )sgn(Y i -Y„) 
n ^ 

js E (R,„-^)(S (i| -2ji); 

11 1 

^TT) Esgn(Xi-Xj) sgn(Yi-Yj); 
V i # j 


and 


{ °> 


n (n — l)(n — 2 ) . ^ . S S n ( X i ~ X j) S g n (Yi-Y k ), fom>2 

j^ k 

» T k 

for n=2, 


r - 1 if u < 0 

where R,^ is the rank of X i5 S/;\ is the rank of Y ; , and sgn(u) = \ 0 if u=0 

1 ' ’ L +1 if u > 0 


We also have 


k' = 


(n-2)k +3t 


The grade correlation coefficient k can be expressed as, see Konijn [9], 

= joo JOO F(x y) dFi(x) dF2(y)) 


(5) 

(6) 


(7) 


( 8 ) 


-oo -OO 


(9) 
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where F(x,y) is the bivariate probability distribution function, and F x (x), F 2 (y) are 
the marginal distributions. 

For the Gumbel a -case, the right side of (9) becomes 
I = J°° exp[-x-y-2 e‘ x -2 e‘ y + a(e x + e y )" 1 ] dx dy 

“OO “OO 

= /°° /°° exp[-x-y-2 e' x (l+ ? =4) + a e‘ x (l+^)" 1 ] dx dy. (10) 

-oo “OO e e 


Let e y ’ x = e l or y-x=t, y=t+x, and (10) becomes 
I = J°° e _t f°° exp[-2x-e" x (2(l+e’ t ) -a(l+e t )' 1 )] dx dt. 

-oo -oo 

The inner integral is 

joo e -2x-Ae x w h ere ^ = 2(l+e" t ) -a(l-fe t )” 1 

-OO 


i f< 
A J 


d(e‘ Ae ) 


L [e" x e’ Ae 


+ f°° e‘ x e‘ Ae 


dx] 


= S [0 + X /°° d( e Ae " X )] = 72 • 

-oo A 

Therefore, 

I — f°° * 1 rri-o dt 

-oo {2(l+e *)-*(!+«) l } 2 

= r°° (1+e \ 2 dt 

-oo {2(l-!-e' t )(l+e t )-a} 2 
— t / ^ . t\2 

= J°° 6 t t 2 dt ‘ 

-oo {2(2+e t +e ^-a} 2 


Let l+e t =2, e fc dt=dz, then 


I - / 


-2t 2 j 

e z dz 


1 {2(2 e l +e zt +l e t )/e t > 2 


= / 


1 {2(l+e l ) 2 -f( Z -l)} 2 

4 f i 2 <J a, 2 
1 - 2 z + 2 > 


1 r oo z^ dz 



Thus equation (9) can be written as 


1 +5= /°° 
1 


z z dz 

(z 2 -fz + f ) 2 


:= I . 


( 11 ) 


^an if q > 0, 
l, if q— o, 


where q=a(2 -|), 0 < a < 1. 
Therefore, 


1+1 = 


\ 1 +>jq^ — tan ^q), ifq>0, 

[1, otherwise. 


( 12 ) 


A table is given in Gumbel and Mustafi [5], which shows the grade correlation 
coefficient k as a function of a (k is almost a linear function of a). Hence a can be 
estimated from the sample grade correlation coefficient using (12). 


B.2. (n) Estimating the parameter m using the standard deviation of the difference 
of the reduced variates 

The quadrant method was based on the value of the critical content Q x 3 . The 
conditions 0 < a,^ < 1 on the parameters impose conditions on Q x 3 , and hence some 
restrictions on the quadrant method. Another method is available to estimate the 
parameter m. 

The difference of the two reduced variates for the bivariate extreme value 
distribution is a logistic distribution (this does not hold for the two dependent reduced 
variates in the a -case). 

Let t=Y-X, X and Y are reduced variates, then the probability distribution of t is 
F(t) = (1 + e mt y\ 

The t variate has zero mean and standard deviation 
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This leads to the following estimate of m 


m = 


^3S t ’ 


(14) 


where S t is the sample value of <r t (obtained from the observed reduced differences). 

Also, it is known that the correlation coefficient p between the reduced variates 
X and Y is related to m by 

m = (l-pp. (15) 

This gives another estimate of m 


m 


(i-p)‘ 


(16) 


B.3. The “strip estimate” 

Posner et al. [11] proposed a method of estimating a or m based on the 
distribution of |X-Y|. Recall that in general the bivariate extreme value distribution, 
see equation (III. 3), is given by 

F(x,y) = exp[- (e" x + e' y ) g(y-x)], 
where g is the dependence function. 

In this case they proved that 

P {a<X-Y<b} = h(b) - h(a) , (17) 

where 

h(t) = (1 + e*)-' + ^ (18) 

In particular, if g is given by 
g(t) = 1 - f sech 2 (|) := 1 - | u(t) , 

then 

P{|X-Y|<<) = + := P(<S) . 

If is the number of times |X — Y| < 6 in a sample of size N, then is 

an estimator of the probability P(£) above. Thus an estimator a of a is obtained by 
solving the equation 
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i'r#) _ e^-l , 2g'((5) 

N “ e«+l *C«) 

= e^. 2 <l> "'<*> 

J+l ! + (f) “(*) ’ 

where 


u(t) = sech 2 (I), and 
u'(t) = -sech 2 (|) tanh(|). 


For the m-case: 


m t i i \m 


g(t)=^ ) 


g'(t) = 


, 1 < III < OO. 

t i i \/_mt , i \rn _ ^ „mt / mt . 1 \fn _t 


(e l +l)(e mi +l) m e mi - (e mi +l) J 

(«‘+l) 2 


g'(t) _ e mt 
g(t) 


mt 


mt . 1 mt . ^ 

e +1 e +1 


Therefore, an estimator m of a is obtained by solving the equation 
^(6) 


e - 1 

6 , i 

e +1 


_ / rii<5 £> \ 

+ 2 fer: "fr) 


'e ri,,5 +t e 5 +l' 


n m 6 

2e -l 


e^+l 


e™ 6 -l 

e^+l 


(19) 


(20) 
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VI. HYPOTHESES TESTING AND GOODNESS-OF-FIT TEST 

Tests of hypotheses on parameters can be done by relating the acceptance region to 
the confidence region. Confidence interval estimation is available in the literature, see 
for example [19], pp. 49-80. 

A goodness-of-fit test for the smallest extreme is Gumbel distribution against the 
alternative, which is Weibull extreme distribution; this was discussed by Montfort [10]. 

Another test of the hypotheses that the largest extreme is Gumbel against the 
alternative that it is either Weibull or Frechet is due to Tiago de Oliveira [17,18]. 

As pointed out in B.2.(ii), if t=Y-Y, where X and Y are dependent reduced 
variates of the Gumbel type, then t has a logistic distribution in the m - case but not 
the a -case. This fact could be used to discriminate between the choice of the a -case 
and the m - case. 
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